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ACCELERATION OF COSMIC RAYS AT COSMIC SHOCKS 
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ABSTRACT 

Nonthermal particles can be produced due to incomplete thermalization at collisionless shocks and further 
accelerated to very high energies via diffusive shock acceleration. In a previous study we explored the cosmic ray 
(CR) acceleration at cosmic shocks through numerical simulations of CR modified, quasi-parallel shocks in ID 
plane-parallel geometry with the physical parameters relevant for the shocks emerging in the large scale structure 
formation of the universe (Kang & Jones 2002). Specifically we considered pancake shocks driven by accretion 
flows with u Q = 1500 km s _1 and the preshock gas temperature of To = 10 4 — 10 8 K. In order to consider the 
CR acceleration at shocks with a broader range of physical properties, in this contribution we present additional 
simulations with accretion flows with u = 75 — 1500 km s _1 and To = 10 4 K. We also compare the new simulation 
results with those reported in the previous study. For a given Mach number, shocks with higher speeds accelerate 
CRs faster with a greater number of particles, since the acceleration time scale is t acc oc u~ 2 . However, two 
shocks with a same Mach number but with different shock speeds evolve qualitatively similarly when the results 
are presented in terms of diffusion length and time scales. Therefore, the time asymptotic value for the fraction 
of shock kinetic energy transferred to CRs is mainly controlled by shock Mach number rather than shock speed. 
Although the CR acceleration efficiency depends weakly on a well-constrained injection parameter, e, and on shock 
speed for low shock Mach numbers, the dependence disappears for high shock Mach numbers. We present the 
"CR energy ratio", <&(M S ), for a wide range of shock parameters and for e = 0.2 — 0.3 at terminal time of our 
simulations. We suggest that these values can be considered as time-asymptotic values for the CR acceleration 
efficiency, since the time-dependent evolution of CR modified shocks has become approximately self-similar before 
the terminal time. 
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I. INTRODUCTION 

Collisionless shocks form ubiquitously in astrophys- 
ical environments via collective electromagnetic vis- 
cosities, when supersonic disturbances propagate into 
tenuous, magnetized, cosmic plasmas. Due to in- 
complete plasma thermalization at collisionless shocks, 
some suprathermal particles leak upstream and their 
streaming motions against the background fluid gener- 
ate strong MHD Alfven waves upstream of the shock 
(e.g., Wentzel 1974; Bell 1978;Lucek & Bell 2000). Al- 
though these self-excited MHD waves provide necessary 
electromagnetic viscosities to confine thermal particles 
to the downstream region of the shock, some suprather- 
mal particles in the high energy tail of the Maxwellian 
velocity distribution may re-cross the shock upstream. 
Then these particles are scattered back downstream by 
those same waves and can be accelerated further to 
higher energies via Fermi first order process. Here we 
refer to cosmic rays as nonthermal particles above the 
Maxwellian distribution in momentum space, so they 
include nonrelativistic, suprathermal particles as well 
as usual relativistic particles. It is believed that cosmic 
ray particles are natural byproducts of the collisionless 
shock formation process, and they are extracted from 
the shock-heated thermal particle distribution (Malkov 
& Volk 1998, Malkov & Drury 2001). 



According to diffusive shock acceleration (DSA) the- 
ory, a significant fraction of the kinetic energy of the 
bulk flow associated with a strong shock can be con- 
verted into CR protons, depending on the CR injec- 
tion rate (e.g., Drury 1983; Berezhko, Ksenofontov, & 
Yelshi 1995; Kang & Jones 1995). In Gieseler, Jones 
& Kang (2000), we developed a numerical scheme that 
self-consistently incorporates a "thermal leakage" in- 
jection model based on the analytic, nonlinear calcula- 
tions of Malkov (1998). This injection scheme was then 
implemented into the combined gas dynamics and the 
CR diffusion-convection code with the Adaptive Mesh 
Refinement technique by Kang, Jones & Giesler (2002). 
In Kang & Jones (2002) (Paper I, hereafter) we applied 
this code to study the cosmic ray acceleration at shocks 
by numerical simulations of CR modified, quasi-parallel 
shocks in ID plane-parallel geometry with the physi- 
cal parameters relevant for the cosmic shocks emerging 
in the large scale structure formation of the universe. 
We adopted the Bohm diffusion model for CRs, based 
on the hypothesis that strong Alfven waves are self- 
generated by streaming CRs. We found in Paper I that 
about 10~ 3 of incoming thermal particles are injected 
into the CR population with our thermal leakage in- 
jection model and up to 60 % of initial shock kinetic 
energy is transferred to CRs for strong shocks. 
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During the formation of large scale structure in the 
universe cosmic shocks are produced by flow motions 
associated with the gravitational collapse of nonlinear 
structures (e.g., Kang et al. 1994a; Miniati et al. 2000). 
These structures are surrounded by accretion shocks 
and CRs can be accelerated to very high energies at 
these shocks via the first order Fermi process (Kang, 
Ryu & Jones 1996, Kang, Rachcn, & Biermann 1997, 
Miniati et al. 2001a, b). In a recent study, Ryu, 
Kang, Hallman, & Jones (2003) studied the charac- 
teristics of cosmic shocks in a cosmological hydrody- 
namic simulation of a ACDM universe and their roles 
on thermalization of gas and acceleration of nonther- 
mal particles. They showed that cosmic shocks have a 
wide range of physical parameters with shock speed 
up to ~ 3000 km s , the preshock temperature of 
10 4 < T < 10 8 , and Mach number up to a few 100. 
In this contribution we extend our study in Paper I by 
including simulations of shock models with a broader 
range of physical parameters that are relevant for such 
cosmic shocks. 

In the next section we briefly describe our numerical 
simulations. The simulation results are presented and 
discussed in §111, followed by a brief summary in §IV. 

II. Numerical Simulations 
(a) Numerical Methods 

Our numerical code was described in detail in Paper 
I, so here we briefly summarize some special features 
that are different from standard hydrodynamics codes: 
1) We solve the gasdynamic equations with CR pres- 
sure terms added in the conservative, Eulerian formula- 
tion for one dimensional plane-parallel geometry along 
with the diffusion-convection equation for the CR mo- 
mentum distribution function (Kang et al. 2001). 2) 
An adaptive mesh refinement technique is applied to 
shock fronts that are tracked as discontinuous jumps by 
a front-tracking method (Kang et al. 2001). 3) An ad- 
ditional equation for the "Modified Entropy" is solved 
to follow accurately the adiabatic changes outside of 
shocks, particularly in the precursor region (Kang et 
al. 2002). 4) We adopt an injection scheme based on 
a "thermal leakage" model that transfers a small pro- 
portion of the thermal proton flux through the shock 
into low energy CRs (Gieseler et al. 2000). There 
is a free parameter in the adopted injection model: 
e = Bq/B±, defined to measure the ratio of the am- 
plitude of the postshock MHD wave turbulence B± to 
the general magnetic field aligned with the shock nor- 
mal, Bo (Malkov and Volk, 1998). This parameter is 
rather well constrained, since 0.3 <, e < 0.4 is indicated 
for strong shocks (Malkov & Volk 1998). However, such 
values lead to very efficient initial injection and most of 
the shock energy is transfered to the CR component for 
strong shocks of M s ^ 30 (Kang et al. 2002), causing 
a numerical problem at the very early stage of simula- 
tions. So in this study we consider models with e = 0.2 



and 0.3 and then discuss how the CR acceleration de- 
pends on e. 

(b) Diffusion Coefficient, Time, and Length 
Scales 

As in Paper I we assume a Bohm type diffusion 
model in which the diffusion coefficient is given as 

P 2 

<P) = K ° {p 2 + 1 y/2 > (!) 

where p is expressed in units of m p c, n Q = 3.13 x 
10 22 cm 2 s _1 i?~ 1 , and B^ is the magnetic field strength 
in units of microgauss. This assumption is based 
on the hypothesis that strong Alfven waves are self- 
generated by streaming CRs and provide random scat- 
terings strong enough to scatter particles within one 
gyration radius. In order to model amplification of 
self-generated turbulent waves due to compression of 
the perpendicular component of the magnetic field, 
the spatial dependence of the diffusion is modeled as 
n(x,p) — k(p)po/p(x), where po is the upstream gas 
density. 

With the Bohm diffusion model the mean accelera- 
tion time scale for a particle to reach momentum p is 
related to the diffusion coefficient and shock speed as 

race « ^ « (7.9 x 10 9 years)(p 10 )S;V 1 o 2 , (2) 

* 8 

where p w = p/10 10 and y 100 o = V^/IOOO km s" 1 is 
shock speed in units of 1000 km s _1 . Although this 
expression is valid for strong shocks, it gives a reason- 
able estimate for any shock strength within a factor of 
a few (see Paper I). The so-called diffusion length of 
CR protons is given by 

Aiiff = ^ « (l.OMpc^oB^ooV (3) 

Thus, for a given value of n , the CR acceleration pro- 
ceeds faster and CR particles diffuse on smaller length 
scales at shocks with greater speeds. Consequently, it 
is more convenient to normalize the simulation vari- 
ables in terms of a diffusion time, t = k /u^, and a 
diffusion length, x = k /u (where u is a charac- 
teristic velocity in the problem), rather than adopting 
constant values across different models. They represent 
the diffusion length and time scales for the protons of 
p = 1. In our discussion a choice of n (or B^) is arbi- 
trary, since we focus on the CR acceleration efficiency 
at time asymptotic limits. 

(c) ID Plane-parallel Shock Models 

In general, cosmic shocks associated with the large 
scale structure formation can be oblique and have 
various geometries, depending on types of nonlinear 
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structures onto which accretion flows fall. As in Pa- 
per I, however, we study the CR acceleration at one- 
dimensional (ID) quasi-parallel shocks which form by 
accretion flows in a plane-parallel geometry in this pa- 
per. Since our simulations follow the acceleration of 
CRs up to only mildly relativistic energies, the diffusion 
length scales are much smaller than typical curvature 
radius of cosmic shocks. So, the diffusion and accel- 
eration of CRs can be studied with fD plane-parallel 
shock models. 

We consider a pancake shock formed by a steady ac- 
cretion flow with a constant density and gas pressure: 
a ID simulation box with [0, x max ] and an accretion 
flow entering into the right boundary of the simulation 
box with a constant density, po, gas pressure, T gl o, and 
velocity, u - There are no pre-existing CRs in the ac- 
cretion flows. The accretion flow Mach number is de- 
fined by M = |?io|/c s , where c s = ("fPgfi/po) 1 ^ 2 is the 
sound speed of the accreting gas. The flow is reflected 
at the left boundary (x — 0) and a shock propagates to 
the right. For a hydrodynamic shock without the CR 
pressure, the shock speed is u s = \uo\/(r— 1) in the sim- 
ulation frame and u' s = \uo\r/(r— 1) in the far upstream 
rest frame, where r is the compression ratio across the 
shock. For the accretion Mach number, 2 < Mq < 100, 
the initial shock speed ranges u' s q — (4/3- 3/2) Kl 
before the CR pressure is built up. So the initial shock 
Mach number is given by M s = (4/3 - 3/2)M . After 
we include the CR acceleration, the CR pressure feed- 
back slows down the shock, so the instantaneous shock 
speed, u' s (t), in the far upstream rest frame decreases 
over time. A CR modified shock consists of a smooth 
precursor and a subshock, since CRs diffuse upstream 
of the subshock, and the CR pressure decelerates and 
heats the preshock flow adiabatically, resulting in weak- 
ening of the subshock. 

(d) Constant u versus Constant T Models 

In Paper 1 we considered a set of models in which 
the infall speed was fixed at uq = —1500 km s~ 4 , 
while the preshock temperature was varied according to 
T = f 4 A(M /100)~ 2 , where the accretion flow Mach 
numbers in the range of 2 < Mo < 100 were consid- 
ered. So, in Paper I we calculated CR modified shocks 
with different preshock temperature but with the same 
accretion speed. We will call this group of models "con- 
stant uq" models, which was designed to study shocks 
with a wide range of the preshock temperature and a 
canonical velocity inside intracluster medium of galaxy 
clusters. 

In this paper we consider accretion flow models with 
a same preshock temperature at To = 10 4 but with 
different shock speeds according to 

K| = 1500 km s- 1 (Mo/100), (4) 

which will be referred as "constant To" models. This set 
of the constant To model is designed to understand the 
CR acceleration at shocks that have a same preshock 



temperature but a wide range of shock velocity. The 
preshock temperature To = 10 4 is chosen to represent 
the mean temperature of the photoionized plasma after 
reionization of intergalactic medium at z ~ 5 (Gncdin 
2000). We consider models with 5 < M < 50, cor- 
responding to 75 km s _1 < |ito| < 750 km s _1 . The 
model with Mo = 100 of the constant To models is 
identical to that of the constant u models presented 
in Paper I, so not considered here. 

(e) Normalization of Physical Variables 

Ideal gasdynamic equations in ID planar geometry 
do not contain any intrinsic time and length scales, 
but in CR modified shocks the CR acceleration and 
the precursor growth develops over the diffusion time 
on diffusion length scales. So, there are intrinsic simi- 
larities in the dynamic evolution and structure of two 
CR shock models with a same Mach number but with 
different preshock temperature or with different shock 
speeds, if we present the results in the coordinate sys- 
tem normalized with diffusion scales (t and x ). 

As in Paper I we set the far upstream density and 
flow values as po = po/p = 1, uq = uq/u q = — 1 in 
code units for all models, where p Q and u Q arc nor- 
malization constants for the gas density and the flow 
speed, respectively. The normalization constants de- 
pend on the accretion Mach number Mo as follows: 
u a = 1500 km s -1 (Mo/100) and so /3 = u Q /c = 

0. 005(M /100). Thus the length and time scales de- 
pend on Mo according to x = n /u oc Mq 1 and 
t Q = K /u^ oc Mq 2 , respectively, with an arbitrary 
choice of diffusion coefficient, n a . The gas density 
normalization constant, p , is arbitrary as well, but 
the pressure normalization constant depends on M as 
T = PoU 2 , oc Mq. Throughout the paper and in the 
code physical variables are given in units of the nor- 
malization constants, x Q , t a , u Q , p a , and P Q . 

The models with smaller Mo have smaller shock 
speeds and lower postshock temperature, so typical 
shock thickness, which is of order of gyro radius of 
the thermal protons, is smaller than that of the mod- 
els with larger M . Since an effective numerical shock 
thickness is the grid spacing at the finest grid level, 
for low Mo models, the required grid spacing in cur- 
rent simulations is much smaller than that of Paper 

1. The new simulations are carried out on a base grid 
with Aa;o = 0.002 using l max = 7 additional grid lev- 
els, so AX7 = 1.56 x 10~ 5 at the finest grid level. This 
leads to a severe requirement of computation time, so 
we run the simulations for t/t a — 20, much shorter 
time than that of Paper I, t/t Q = 200. The simu- 
lated space is x = [0, 20] and N = 10000 zones arc 
used on the base grid. The number of refined zones 
around the shock is N r f = 50 on the base grid and 
so there are 2N r f = 100 zones on each refined level. 
To avoid technical difficulties, the multi-level grids are 
used only after the shock propagates away from the left 
boundary at the distance of x s — 0.05. After the shock 
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Fig. 1. — Time evolution of the shocks driven by the accretion flows with Mo = 5. Heavy lines represent the model with 
T = f0 4 K and u a = 75 km s _1 at t/t„ = 5, 10, and 20. Light lines show the model from Paper I with T = 4 x 10 6 K 
and Wo = 1500 km s _1 at t/t = 5, 10, 20, 100, and 200. The diffusion time scale, to = k /Wo> and the diffusion length, 
x = n /u are defined by the accretion speed of each model, u . The inverse wave-amplitude parameter e = 0.2 is adopted 
for both models. The accretion flows are reflected at x/x = and the shocks of M s — 6.8 propagate to the right, so the 
leftmost profile corresponds to the earliest time. Note the distance from the reflecting plane is in a logarithmic scale. 



moves to x s = 0.05 (at t « 0.15 for strong shocks), 
the AMR technique is turned on and the CR injec- 
tion and acceleration are activated. This initial delay 
of the CR injection and acceleration should not affect 
the final outcomes. For all models we use 230 uni- 
formly spaced logarithmic momentum zones in the in- 
terval \og(p/m p c) — [logpoi logpi] = [—3.0, +3.0] As 
in our previous studies, the function g(p) = p 4 f(p) is 
evolved instead of f(p) and y = hip is used instead of 
the momentum variable, p for that step. 

III. RESULTS 

Since the CR acceleration is mainly determined by 
velocity jump across the shock (i.e., Ap/p oc Au/c), 
which is a function of shock Mach number only, we 
expect two models with the same Mach number but 
with different u may have qualitatively similar results, 
when expressed in terms of diffusion scales. The main 



difference between the two models is the "effective " in- 
jection momentum, pinj/rripC ~ 0.01(u o /1500 km s _1 ). 
Here the "effective " injection momentum refers to a 
mean value of the momentum range over which thermal 
leakage takes place (Gieseler et al. 2000). The diffusion 
coefficient of injected, nonrelativistic CRs, n(p) oc pf n ^ 
iov p in j <C 1, has different momentum dependence from 
that of relativistic CRs, n(p) oc p for p > 1. So the sim- 
ilarity between the two models could be broken at early 
evolutionary stages when the CR pressure is still domi- 
nated by nonrelativistic particles. We attempt to make 
comparisons of such two models in this section. 

(a) Modified Shock Structure 

We show the time evolution of the shock driven by 
the accretion flow with M = 5 and u = 75 km s _1 at 
t/t Q = 5, 10, 20 in Figure 1 (heavy lines). The model 
with Mo = 5 and u a = 1500 km s _1 from Paper I is also 
shown at t/t Q = 5, 10, 20, 100 and 200 (light lines). As 
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Fig. 2. — Same as Figure 1 except Mo = 30. Heavy lines represent the model with T = 10 4 K and u a = 450 km s 1 at 
t/t = 5, 10, and 20. Light lines represent the model with T = 4.4 x 10 5 K and u = 1500 km s _1 at t/t a = 5, 10, 20, 100, 
and 200. 



discussed in the previous section, the results are pre- 
sented in two different sets of physical coordinates, so, 
for example, the ratio of diffusion length scales between 
two model is x (75 km s _1 )/a; o (1500 km s _1 ) = 20. 
Since we integrate the model with u Q = 75 km s _1 
up to t /t = 20 when nonrclativistic particles still 
contribute significantly to the CR pressure, we expect 
to see some differences in the two models before that 
time. First of all, in the model with u Q = 75 km s _1 , 
CR particles are injected at much lower injection mo- 
menta and their acceleration time scales in units of 
t a are shorter. So the CR pressure increases faster 
and the modification to the flow structure due the CR 
pressure occurs earlier, compared to the model with 
u = 1500 km s _1 . In both models a substantial pre- 
cursor develops upstream of the shock, resulting in a 
significant modification of flow structure. Also, the 
shock structures seem to have reached approximate 
time-asymptotic states by the terminal time of both 
simulations (t/t a = 20 for u = 75 km s" 1 model and 
t/t Q = 200 for u = 1500 km s" 1 model). Once the 
postshock CR pressure becomes constant, the shock 



structure evolves approximately in a "self-similar" way, 
because the scale length of shock broadening increases 
linearly with time. By comparing the shock struc- 
tures at the terminal times, we expect that the overall 
evolution and shock structure would approach to time 
asymptotic states that are similar for the two models. 

In Figure 2 we make the similar comparison be- 
tween the two models with u = 450 km s _1 and 
u Q = 1500 km s _1 for strong shocks driven by the ac- 
cretion flow with Mq = 30. Due to stronger nonlinear 
feedback effects, the overall shock structure displays 
greater differences at a given value of t/t a , compared 
to the models with M = 5. For example, distances of 
the shock from the reflecting plane are different and the 
gas and CR pressure have very different profiles before 
t/t a < 20. However, the time asymptotic values of the 
CR pressure become similar at the terminal times of 
the two models. 

Figure 3 shows the shock structure of the constant To 
models with different M at the terminal time (t/t a = 
20). Here the spatial coordinate is shifted so the shock 
is located at x/x Q = 0.0 for all models. The degree 
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e = 0.2, t/t = 20., M = 5, 10, 20, 30, 50 




-5 5 10 -5 5 10 




Fig. 3. — Structure of cosmic-ray modified shocks driven by the accretion flows with To = 10 4 K at t/t = 20. The accretion 
flow speed scales with the accretion Mach number as uq — 1500 km s (Mo/100). The inverse wave-amplitude parameter 
e = 0.2 is adopted. The length and time coordinates are expressed in units of x a and to, respectively, defined by the accretion 
speed of each model. 



of nonlinear modification and density enhancement are 
greater in models with higher Mo due to the higher CR 
pressure, so the distance of the shock from the reflecting 
plane (the leftmost point) is smaller. The CR pressure 
seems to converge to a limiting value of P c / p (u' s Q ) 2 ~ 
0.6, where u' s is the initial shock speed before the 
CR modification slows down the shock. The effects of 
nonlinear feedback such as the preshock deceleration 
and the compression through the total shock transition 
can be compared most clearly in the flow velocity plot. 

(b) Particle Distribution Function 

Figure 4 shows the evolution of the particle distribu- 
tion functions (including both thermal and nonthermal 
populations) at the shock (i.e., f[x s ,p]) for two sets 
of models (M = 5 and M = 30) shown in Figures 
1-2. In addition to the usual form of distribution func- 
tion, g(p) = ,f(p)p 4: , we also plotted f(p)p 3 to show the 



differential number density of CRs in logarithm mo- 
mentum bin, dN = f(p)pdhip. In fact g(p) repre- 
sent the differential energy density of CRs in logarithm 
momentum bin, dE — f(p)p 4 dlnp for relativistic mo- 
menta (p ^> 1). First we note that peak momentum 
of thermal distribution is lower for the lower postshock 
temperature, so the particle injection starts at lower in- 
jection momenta in the models with lower u . For the 
model with M — 30 and u a = 1500 km s _1 , we can 
see that the Maxwellian distribution shifts to lower mo- 
menta as the postshock temperature diminishes due to 
energy transfer to CRs. For other models this adjust- 
ment happens well before t/t a = 5, so the postshock 
thermal distribution is almost steady after that time. 

We identify p max as the momentum above which 
g(p) drops sharply, characterizing the effective upper 
cutoff in the CR distribution. This momentum is ap- 
proximately related to the age of the shock in units of 
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Fig. 4. — Evolution of the particle distribution function at the shock, represented as p 4 f(p) and p 3 f(p), is shown for two 
sets of the models shown in Figures 1-2. Heavy lines represent the shock models with To = 10 4 K at t/t a = 5, 10, and 20. 
The accretion flow speed is u a — 75 km s _1 for Mo = 5 and u a — 450 km s _1 for Mo = 30. Light lines represent the shock 
models with u a — 1500 km s _1 at t/t = 5, 10, 20, 100, and 200. The peaked distribution at low momenta (p < 10~ 2 ) 
represent the thermal Maxwellian distribution. The same line types are used here as in Figures 1-2. 



the diffusion time as p max ~ 4.5(t/t ) for p max > 1. 
This explains why values of p max are similar at a given 
value of t/t for all models, regardless of u , injec- 
tion momentum, or shock Mach number. In com- 
parison of the two models with a same Mo, the one 
with lower u a injects CRs at lower momenta, so, for 
example, logp; n j w —3.3 for u Q = 75 km s _1 , while 
logpi n j w —2 for u Q = 1500 km s _1 . Just above 
the injection pool, the distribution function changes 
smoothly from the thermal distribution to an approxi- 
mate power-law whose index is close to the test-particle 
slope for the subshock, i.e., q s = 3r s /(r s — 1) (where r s 
is the compression ratio across the subshock). The dis- 
tribution function g{p) shows the characteristic "con- 
cave upwards" curves reflecting modified shock struc- 
ture (including the precursor) for all models. Although 
the distribution function at a given momentum is larger 



in the higher u models, but the ratio of P c /p u^ is 
larger in lower u models. This leads to slightly greater 
nonlinear feedback in lower u a model, resulting in more 
concave curves. 

Figure 5 shows the total CR distribution within the 
simulation box, G(p) = f p 4 / cr (p)dx and its power law 
slope q = — (d\nG/d\np — 4) at the terminal time of 
three sets of simulations. The models shown in top 
two panels are the constant To models with e = 0.2 
for Mo = 5, 10, 20, 30, and 50. The particle momenta 
near logp w —3.3 are injection momenta above the 
thermal distribution and have similar values within a 
factor of two for models with different M , indicating 
similar thermal populations. Although the postshock 
temperature increases with shock speed as T 2 oc (u s t) 
in pure gasdynamic shocks, all models in this group 
have similar T 2 . This is because the postshock gas 
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Fig. 5. — CR distribution function integrated over the simulation box, G(p) = J p 4 / cr (p)dx, and its power law slope, 
q = — (d In G/d hip — 4), at the terminal time of each simulation are shown. Top two panels are for the models with To = 10 4 
K and e = 0.2, while middle two panels are for the models with To = 10 4 K and e = 0.3. For both sets of models the 
accretion speed is given as u = 1500 km s _1 (Mo/100) and the terminal time is t/t = 20. Bottom two panels show the 
models with u a — 1500 km s _1 and e = 0.2. The preshock temperature is given by To = 10 4 K(Mo/100) -2 . The curves for 
G(p)are labeled with the accretion Mach number. 



pressure falls and the postshock density increases due 
to nonlinear feedback to a greater extent at higher 
Mach number models. The integrated distributions 
also show the characteristic "concave upwards" curves 
and this "flattening" trend is stronger for higher Mq 
models. The slope of the total CR spectrum ranges 
over 4.0 <> q ^ 4.7 near p m j then decreases with the 
particle momentum and converges for strong shocks 
to q <~ 2 just below p max . The same kind of mod- 
els but with e = 0.3 are shown in middle two panels. 
They show that the injection takes place at slightly 
lower momenta with greater numbers of particles in 
the Maxwellian tail for the models with larger e, lead- 
ing to the slope near pj n j as large as q ~ 5. So the CR 
distribution functions for e = 0.3 models are greater 
than those for e = 0.2 models for p < 1, while they are 



roughly similar for p > 1 . One can expect that the CR 
pressure of these two sets of models may be different 
at early stage when nonrclativistic, fresh injected parti- 
cles dominate the pressure, but it becomes similar once 
relativistic particles dominate. As in Figure 4, all mod- 
els shown here have similar values of p max regardless of 
values of u , since the results are shown at the same val- 
ues of t/t - The bottom two panels show the constant 
u Q models from Paper I. Since u a = 1500 km s _1 for 
all models, the postshock temperature would be simi- 
lar if these are pure gasdynamic shocks. Due to their 
higher degree of nonlinear modification to the struc- 
ture, however, the postshock thermal populations have 
lower temperature in the models with higher Mo. 

Although the evolution of these shock becomes ap- 
proximately self-similar and the postshock CR pressure 
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t/l. t/t. 

Fig. 6. — The ratio of total CR energy in the simulation box to the kinetic energy in the initial shock frame that has 
entered the simulation box from upstream, 3>(i), the postshock CR pressure in units of far upstream ram pressure in the 
instantaneous shock frame, and time-averaged injection efficiency, £(t). Left three panels are for Mo = 5 — 50 and e = 0.2. 
Right three panels show the same quantities for Mo = 5 — 50 and e = 0.3. 



reaches a quasi-steady value before the terminal time, 
CR particles continue to be accelerated to ever higher 
energies and so the CR distribution continues to ex- 
tend to higher momenta. Thus the CR distribution 
functions in Figure 5 show only a snap shot at the ter- 
minal time. But they illustrate how the CR distribu- 
tion deviates from the test-particle like power-law due 
to nonlinear feedback. 

(c) Injection and Acceleration Efficiencies 

As in Paper I we dehne the injection efficiency as 
the fraction of particles that have entered the shock 
from far upstream and then are injected into the CR 
distribution: 



£(*) 



Jo 



dx J^4TrfcR(p,x,t)p 2 dp 



(5) 



where fen is the CR distribution function, no is the 
particle number density far upstream, u' s is the instan- 



taneous shock speed, and t\ is the time when the CR 
injection/acceleration is turned on. 

As a measure of acceleration efficiency, wc define the 
"CR energy ratio" ; namely the ratio of the total CR en- 
ergy within the simulation box to the kinetic energy in 
the initial shock frame that has entered the simulation 
box from far upstream, 



$(t) = ^ 



dx£; CR .(x, t) 



0.5po«o) 3 i 



(6) 



where u' s is the initial shock speed. Since the shock 
slows down due to nonlinear modification, the kinetic 
energy flux in the instantaneous shock rest frame also 
decreases. So the ratio of the total CR energy in the 
simulation box to the kinetic energy defined in the in- 
stantaneous shock frame which has entered the shock 
from far upstream, i.e., 



r 



dxE C R(x,t) 



J 0.5p o u o (i) 3 dt 



(7) 
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Fig. 7. — The CR energy ratio, <J?, at the terminal time of the three sets of simulations shown in Figure 5 as a function 
of the shock Mach number, M s . The solid line with circles are for models with u a = 1500 km s~ 4 and e = 0.2, the dotted 
line with triangles for models with To = 10 4 K and e = 0.2, and the dashed line with squares for models with To = 10 4 K and 
e = 0.3. 



can be much larger than especially for high Mach 
number shocks. 

Figure 6 shows the CR energy ratio, the CR pres- 
sure at the shock normalized to the ramp pressure of 
the upstream flow in the instantaneous shock frame, 
Pc.2 / Po{ u ' s ) 2 , and the "time-averaged" injection effi- 
ciencies, £, for models with different M when e = 0.2 
(left three panels) and e = 0.3 (right three panels). 
For all Mach numbers the postshock P c .2 increases un- 
til a balance between injection/acceleration and advec- 
tion/diffusion of CRs is achieved, and then stays at a 
steady value afterwards. The time-asymptotic value 
of the CR pressure becomes P c ,2/po( u ' s o) 2 ~ 0.6 in 
the initial shock frame and P c ,2 / Po(u' s ) 2 ~ 0.9 in the 
instantaneous shock frame for Mq = 50 with e = 0.2. 
After P c _-i has reached a quasi-steady value and the evo- 
lution of the P c spatial distribution has become "self- 
similar" , the CR energy ratio also asymptotes to a con- 
stant value. Time-asymptotic value of $ increases with 
M , but it converges to $ « 0.53 for M > 20 and 
e = 0.2. As discussed in detail in Kang et al. (2002), the 
injection rate is higher for higher subshock Mach num- 
ber and for larger values of e, which leads to the higher 
CR pressure and higher <J>. Once again, however, this 
dependence on e becomes weaker for stronger shocks. 
The average injection rate is about £ ~ 10 1 — 10 -3 8 
with e = 0.2 and £ w 10~ 3 with e = 0.3. 

Figure 7 shows values of the CR energy ratio $ at 
the terminal time of three sets of simulations shown in 
Figure 5. Although we were able to follow CR accelera- 
tion only up to p m ax ~ 10 — 100, these values can serve 
as estimates for the time-asymptotic CR acceleration 
efficiency, since $ seems to approach constant values 
after t/t a > 1. As shown in Figure 6 the CR acceler- 
ation is more efficient for larger e, so, for the constant 



T models with M s = 6.8, $ « 0.39 for e = 0.2 (dot- 
ted line), while $ « 0.44 for e = 0.3 (dashed line). 
The difference becomes smaller for high Mach num- 
ber shocks (Ms <^ 30). The constant T models have 
smaller u Q = 1500 km s _1 (100/M ) for M < 100 than 
the constant u D = 1500 km s _1 models. Due to this ve- 
locity difference, for given values of Mach number and 
e, the constant T models have slightly higher $ than 
the constant wo models. For the constant uq model 
with M s = 6.8, for example, $ w 0.32 for e = 0.2 (solid 
line). 

IV. SUMMARY 

We have calculated the CR acceleration at ID quasi- 
parallel shocks by using our cosmic-ray AMR Shock 
code (Kang et al. 2002), which incorporates the "ther- 
mal leakage" injection process to the CR/hydro code 
that solves the CR diffusion-convection equation along 
with CR modified gasdynamic equations. Our sim- 
ulations have been performed in a ID plane-parallel 
space in which shocks are driven by the accretion flows 
with u Q = 75 — 1500kms _1 and the temperature of 
T = 10 4 K. Mach number of the resulting shocks ranges 
6.8 < Ms < 133. We have compared the new simula- 
tion results with those from the previous study in which 
the accretion flows with u Q = 1500kms _1 and the tem- 
perature of T = 10 4 — 10 8 K were considered (Kang & 
Jones 2002, Paper I). 

Detailed simulation results found in Paper I remain 
valid, so we briefly review the main conclusions here to 
make the present paper self-contained. 

1) Suprathermal particles can be injected very effi- 
ciently into the CR population via the thermal leakage 
process, so that typically a fraction of 10~ 4 — 10~ 3 of 
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the particles passed through the shock become CRs for 
e = 0.2-0.3. 

2) For a given value of shock Mach number, the in- 
jection efficiency and the CR acceleration are higher 
for larger e. But this dependency is weaker for higher 
Mach numbers of M s £ 30. 

3) For a given value of e, the acceleration efficiency 
increases with M a , but it asymptotes to a limiting value 
of the CR energy ratio, $ w 0.5, for M s J> 30 and 
e = 0.2 - 0.3. 

4) The CR pressure seems to approach a steady- 
state value in a time scale comparable to the accelera- 
tion time scales for the mildly relativistic protons after 
which the evolution of CR modified shocks becomes ap- 
proximately "self-similar". This feature enables us to 
predict time asymptotic values of the CR acceleration 
efficiency through numerical simulations of CR shock 
models with a broad range of the physical parameters 
in this work. 

5) We suggested that the CR acceleration is innate 
to collisionlcss shock formation process and CRs can 
absorb a significant fraction of initial shock kinetic en- 
ergy 

The main purpose of this comparison study is to ex- 
plore how the CR acceleration depends on the preshock 
temperature and shock speed as well as shock Mach 
number. Unlike pure gasdynamic shocks, CR modified 
shocks have intrinsic scales, that is, a diffusion time, 
t = Ko/"oi an d a diffusion length, x Q = n /u , because 
the mean acceleration time scale is t acc oc k(p)/u 2 . So 
the CR acceleration and the evolution of CR modi- 
fied shocks proceed faster for models with higher shock 
speeds. On the other hand, two models with a same 
Mach number but with different accretion speeds should 
show qualitatively similar results, when the dynamical 
evolution is presented in terms of diffusion scales, t a 
and x , since the CR acceleration is controlled mainly 
by the velocity jump across the shock transition, which 
depends on shock Mach number. 

For models with lower shock speed, the mean veloc- 
ity of thermal particles is smaller and so the effective 
injection momentum, p; n j, is smaller, compared with 
those for higher shock speed models. Nonrelativistic 
suprathermal particles (pinj < p < 1) are accelerated 
quickly to relativistic energies because of small diffusion 
coefficients (k(p) oc p 2 ). On the other hand, the high- 
est momentum of the CR distribution increases with 
the age of the shock as p max ~ 4.5(i/t ), and so the 
CR distribution function extends to a similar value of 
Pmax at a given t/t for all models. As a result, for 
the model with lower shock speed, the ratio of P c / p a u 2 
increases faster, nonlinear modification to the under- 
lying flow sets in earlier, and the CR energy ratio is 
larger. However, we expect these differences would be- 
come smaller for t/t a 3> 1 when relativistic particles 
dominate the CR pressure. For example, for the two 
models with the same Mach number at M s = 6.8 and 
the same injection parameter at e = 0.2, <f> = 0.32 for 



u D = 1500 km s _1 and $ = 0.39 for u D = 75 km s _1 at 
the terminal time. For larger values of e, $ arc higher 
due to the higher injection rate. For example, $ = 0.44 
for the model with M s = 6.8, u = 75 km s _1 , and 
e = 0.3. Thus, for models with M s = 6.8, the CR en- 
ergy ratio ranges 0.32 < $ < 0.44, depending on the 
shock speed and the injection parameter. The differ- 
ence becomes smaller at higher shock Mach numbers. 
The CR energy ratio at the terminal time of our sim- 
ulations is summarized for different models in Figure 
7. 
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